function actionFrequency()
% =========================================================================
% Shows the problems with the frequency of action
% =========================================================================

% Load file

edff  = ASCFILE('sample.asc');
analysis = ET_ANALYSIS(edff);
analysis.refresh_rate = analysis.resample;


[vd,s,ox,ow] = analysis.normalizedData();

% Transformed to degrees

v2d = @(x) sign(x).*acosd(16./(4 * sqrt(16 + x.^2)));

vd = v2d(vd);
s  = v2d(s);
ox = v2d(ox);
ow = v2d(ow);

nc = 15;

%nc = 1;

vd = vd/nc;
s = s/nc;
ox = ox/nc;
ow = ow/nc;
c = -8/nc:1/nc:8/nc;
%c = -8:8;
% INFERENCE MODEL: DEM.M and DEM.G
%==========================================================================

% create illustrative (M,G) if not specified (EXAMPLE)
%----------------------------------------------------------------------

% hidden causes and states
%======================================================================
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor velocity
%   x.x(1) - target angle - extrinsic coordinates
%
% v    - causal states: force on target
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception)
%   g(2) - oculomotor velocity
%   g(:) - visual input - intrinsic coordinates
%----------------------------------------------------------------------

%%
% Set-up
%======================================================================
M(1).E.s = 1/2;                               % smoothness
M(1).E.n = 4;                                 % order of
M(1).E.d = 1;                                 % generalised motion


% sensory mappings with and without occlusion
%----------------------------------------------------------------------


% Generative model (M) (sinusoidal movement)
%======================================================================
% Endow the model with internal dynamics (a simple oscillator) so that is
% recognises and remembers the trajectory to anticipate jumps in rectified
% sinusoidal motion.

% slow pursuit following with (second order) generative model
%----------------------------------------------------------------------
x.o = [vd(1);0];                                  % motor angle & velocity
x.x = vd(1);                                      % target location


% level 1: Displacement dynamics and mapping to sensory/proprioception
%----------------------------------------------------------------------



%tmpdf = @(x) -2*(c' - (x.x - x.o(1)));
elasticity = 1/20;

%M(1).f = @(x,v,p) [x.o(2); (v - x.o(1))/4 - x.o(2)/2; v - x.x];

%M(1).f = @(x,v,p) [x.o(2); (v - x.o(1))/4 - x.o(2)/2; v ];
M(1).f = @(x,v,p) [x.o(2); v - x.o(2)/2 - x.o(1)*elasticity ; v ];

M(1).g  = @(x,v,p) [x.o(1);x.o(2); exp(-(c' - (x.x - x.o(1))).^2).*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))];

%M(1).g  = @(x,v,p) [x.o(1);x.o(2); exp(-(c' - (x.x - x.o(1))).^2)];

%M(1).g  = @(x,v,p) [x.o(1);x.o(2); tmpf(x)];

%M(1).gx = @(x,v,p) [1,0,tmpf(x)'.*tmpdf(x)'.*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))
%    0,1,zeros(1,17)
%    0,0,-tmpf(x)'.*tmpdf(x)'.*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))]';

%M(1).fx = @(x,v,p) [0,-0.25,0; 1,-0.5,0;0,0,-1]';
%M(1).fv = @(x,v,p) [0;0.25;1];


M(1).x = x;                                   % hidden states
M(1).V = exp(4);                              % error precision observables
M(1).W = diag([exp(4),exp(4),exp(4)]);                              % error precision hidden states


% level 2: With hidden (memory) states
%----------------------------------------------------------------------

M(2).f  = @(x,v,p) [x(2); -x(1)]*v;

M(2).fx = @(x,v,p) [0,1;-1,0]*v;
M(2).fv = @(x,v,p) [x(2); -x(1)];

M(2).g  = @(x,v,P) x(2);
%M(2).gx = @(x,v,P) sparse([1,0]);
%M(2).gv = @(x,v,P) sparse(1,1);

M(2).x  = [0; 0];                             % hidden states
M(2).V  = exp(17)*100;                            % error precision
M(2).W  = exp(4);                             % error precision

% level 3: Encoding frequency of memory states (U)
%----------------------------------------------------------------------
M(3).v = 0;
M(3).V = exp(9.201);


% generative model (G)
%======================================================================

% first level
%----------------------------------------------------------------------

G(1).f = @(x,v,a,p) [x.o(2), a - x.o(2)/2 - x.o(1)*elasticity, x.o(2) ]';
%G(1).fx = @(x,v,a,p) [0,0,0;1,-1/2,0;0,0,-1]';

%G(1).fv = @(x,v,a,p) [0, 0, 1]';
%G(1).fa = @(x,v,a,p) [0, 1, 0]';

%G(1).g = @(x,v,a,p) M(1).g(x,v,p);%@(x,v,a,p) [x.o; exp(-((-8:8)' - ( v - x.o(1))).^2).*~(v > (ox-abs(ow)/2) && v < (ox+abs(ow)/2))];
G(1).g  = @(x,v,a,p) [x.o(1);x.o(2); exp(-(c' - (v - x.o(1))).^2).*~(v > (ox-abs(ow)/2) && v < (ox+abs(ow)/2))];
%G(1).g  = @(x,v,a,p) [x.o(1);x.o(2); exp(-(c' - (v - x.o(1))).^2)];
%G(1).gx = @(x,v,a,p) M(1).gx(x,v,p);
%G(1).ga = @(x,v,a,p) sparse(19,1);

G(1).x = x;                                  % hidden states
G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(4);  % motor gain

% second level
%----------------------------------------------------------------------
G(2).v = 0;                                  % exogenous force
G(2).a = 0;                                  % action force
 

%======================================================================

% Check generative model
%--------------------------------------------------------------------------
M      = spm_DEM_M_set(M);
 

% Experimental input (C)
%==========================================================================

xU.dt = 16.6666; % ms
N     = numel(vd);
xU.u  = vd;


% Priors (U)
%==========================================================================


% (EXAMPLE) prior beliefs about target frequency (w)
%----------------------------------------------------------------------
U = zeros(1,N) + 2*pi/analysis.FpC(1);

% Data and confounds
%==========================================================================

% Generate simulated data (EXAMPLE)
%----------------------------------------------------------------------

% meta-model and true parameters
%----------------------------------------------------------------------
MM.M    = M;
MM.G    = G;
MM.U    = U;


xY.y    = s';

 
% DCT confounds
%==========================================================================
%xY.X0  = spm_dctmtx(N,1);
 
% and [serial] correlations (precision components) AR model
%--------------------------------------------------------------------------
xY.Q   = {spm_Q(1/2,N,1)};
 
 
% META-MODEL parameters
%==========================================================================
% These parameterise the inference scheme function (IS) at the bottom of
% this script - this functions defines how the parameters are used and
% therefore optimised. Change this function to specify the meta-model


% prior expectations (EXAMPLE)
%----------------------------------------------------------------------
%pE.W  = 0.6708;
pE.W  = 10.6;
%pE.W  = -0.6708;
%pE.a  = 0.6137;


% prior covariance (EXAMPLE)
%----------------------------------------------------------------------
pC.W  = 1;
%pC.a  = 1;


% This completes the meta-model specification. The fields are now assembled 
% and passed to spm_nlsi_GN (nonlinear system identification using Gauss-
% Newton-like gradient ascent).
%==========================================================================


% META-MODEL (MM)
%==========================================================================
MM.M   = M;
MM.G   = G;
 
% hyperpriors (assuming about 99% signal to noise)
%--------------------------------------------------------------------------
hE     = 8 - log(var(spm_vec(xY.y)));
hC     = exp(-4);
 
% Meta-model
%--------------------------------------------------------------------------
MM.IS  = @IS;
MM.pE  = pE;
MM.pC  = pC;
MM.hE  = hE;
MM.hC  = hC;

%% Simulations 
%==========================================================================

% By changing the damping, high frequency variations in action are
% prevented making the model more stable.

%damp = 1/5;
MM.pE  = pE;
%G(1).f = @(x,v,a,p) [x.o(2), a - x.o(2)*damp, v - x.x]';
%G(1).fx = @(x,v,a,p) [0,0,0;1,-1*damp,0;0,0,-1]';
%MM.G   = G;
spm_figure('getwin','High Damping');
[Ey,DEM] = IS(pE,MM,xU);
spm_DEM_qU(DEM.qU);


%%
% damp = 1/8;
% 
% G(1).f = @(x,v,a,p) [x.o(2), a - x.o(2)*damp, v - x.x]';
% G(1).fx = @(x,v,a,p) [0,0,0;1,-1*damp,0;0,0,-1]';
% MM.G   = G;
% 
% spm_figure('getwin','Low Damping');
% [Ey,DEM] = IS(pE,MM,xU);  
% spm_DEM_qU(DEM.qU);

keyboard
end

function [y,DEM] = IS(P,M,U)

% parameterise inference model (EXAMPLE)
%--------------------------------------------------------------------------

M.M(2).W = exp(P.W);
%M.G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(P.a);  % motor gain

% reset random number generator (for action schemes)
%--------------------------------------------------------------------------
% Bayesian inference (spm_DEM, spm_LAP, spm_ADEM or spm_ALAP)
%--------------------------------------------------------------------------
DEM.M  = M.M;
DEM.G  = M.G;
DEM.U  = M.U;
DEM.C  = U.u;
DEM    = spm_ADEM(DEM);

% (Hidden) behavioural or electrophysiological response (EXAMPLE)
%--------------------------------------------------------------------------
x      = DEM.Y(1,:);
y      = x';

end